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ABSTRACT 

Context. Using detailed 3D hydrodynamic simulations we study the nature of the Galactic supernova remnant (SNR) CTB 109 
(G109.1-1.0), which is well known for its semicircular shape and a bright diffuse X-ray emission feature inside the SNR. 

Aims. Our model has been designed to explain the observed morphology, with a special emphasis on the bright emission feature 
inside the SNR. Moreover, we determine the age of the remnant and compare our findings with X-ray observations. With CTB 109 we 
test a new method of detailed numerical simulations of diffuse young objects, using realistic initial conditions derived directly from 
observations. 

Methods. We performed numerical 3D simulations with the RAMSES code. The initial density structure has been directly taken from 
^^CO emission data, adding an additional dense cloud, which, when it is shocked, causes the bright emission feature. 

Results. From parameter studies we obtained the position b) - (109.1545°, -1.0078°) for an elliptical cloud with ^cioud = 25 cm"^ 
based on the preshock density from Chandra data and a maximum diameter of 4.54 pc, whose encounter with the supernova (SN) 
shock wave generates the bright X-ray emission inside the SNR. The calculated age of the remnant is about 11,000 yr according to 
our simulations. In addition, we can also determine the most probable site of the SN explosion. 

Conclusions. Hydrodynamic simulations can reproduce the morphology and the observed size of the SNR CTB 109 remarkably well. 
Moreover, the simulations show that it is very plausible that the bright X-ray emission inside the SNR is the result of an elliptical 
dense cloud shocked by the SN explosion wave. We show that numerical simulations using observational data for an initial model can 
produce meaningful results. 
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1. Introduction 

In astronomy, complex morphologies of diffuse objects with fil¬ 
amentary or wispy structures, for example, are often difficult to 
explain because the ambient medium is very inhomogeneous, 
and/or dynamical changes in the flow (e.g. shock waves) modify 
the emission or absorption processes significantly. If one would 
like to know if a certain feature is due to hydrodynamic inter¬ 
actions, a numerical simulation could help to clarify the condi¬ 
tions under which it has originated. A typical case is the X-ray 
observation of supernova remnants (SNRs). Here we present a 
method that exploits the fact that the molecular environment of 
an SNR has a high inertia, and does not evolve much within 
a typical early Sedov-Taylor phase time scale of the order of 
10"^ yr. Therefore, the initial background model, if taken from 
observations, can be scanned identically into a computa¬ 
tional grid. However, the observational background is a projec¬ 
tion onto the plane of the sky and has to be extended to three di¬ 
mensions. Nevertheless, such an initial model represents the real 
situation much better than previous simplifications, such as a ho¬ 
mogeneous background. Density inhomogeneities lead to mass 
loading of the flow and hence to a quite different dynamical and 
thermal history of the plasma. 

The Galactic SNR CTB 109 (G109.1-1.0) associated with 
the anomalous X-ray pulsar (AXP) IE 2259+586 was first de- 
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tected in X-rays by the satellite mission Einstein in 1980 (Gre- 
gory & Eahlman|[l980| ) and has been regularly observed ever 


since. 


Radi o observations at M9 cm ([Hughes et al.||1981[ |1984| ), 
T21 cm ([Hughes et al.||1984[ Braun & Strom||1986|), ^4.6 cm 
( [Hughes et ^.||1984| ), and at 2.7 GHz ( |Downes||1983| ) show a 
morphology that corresponds well to that observed in X-rays 
with the satellite Einstein. However, no radio emission has been 
detected from the AXR The analysis of CO observatio ns (jTatem- 


atsu et al^|1987[ [Kothes et al.||2002] [Taylor et al.||2063| ) has 


revealed that the western part of the SNR has encountered a 
giant molecular cloud (GMC) complex, resulting in a semicir¬ 
cular shape of the remnant. Moreover, in radio a double-shell 
structure can be seen (see Eig.[2(left)). After the first detection 
with Einstein, CTB 109 was observed in X-rays with ROSAT 
PSPC and BBXT (IRho & Petre| 1993 [[ 1997 b, BeppoSAX ([Parmar[ 
eFaLjp^ [1999^ , and A5CA ( [Rho et ai[T99^ . Partial X-ray 


shells were identified in the east, south, and north surrounding 
the AXP, but no shell structure was found in the west owing to 
the GMC. Moreover, an X-ray bright interior region called the 
X-ray lobe was identified in CTB 109. Higher resolution studies 
in X-rays were performed by XMM-Newton ( [Sasaki et al.|2004[), 
Chan dra ( [Sasaki et al.[[2006| [2013[ ), and Suzaku ( [Nakano et~aL] 
2015[ ) to investigate the remnant in more detail. The y-ray ob¬ 


servations by Eermi-hAI ( [Castro et al.[[2012| complement the 
multi-wavelength view of CTB 109. 
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The determination of the distance to an SNR is quite 
challenging. Recent measurements show that the distance to 
CTB 109 l ies between 3 and 4 kpc ( [Kothes et al.||2002 |Tian 


et al. 


|2010| ); the latest value is 3.2 + 6.2 kpc ( [Kothes & 


^oster 


2012). It is even more difficult to determine the proper motion 


of the associated anomalous X-ray pulsar. The analysis by [Ka 
plan et al. ( 2009| ) yields = (-11 ± 20 masyr“^-35 + 

20 masyr“^), which is larger than the v alues (Ury, = (-6.4 + 


0.6 masyr \-2.3 + 0.6 masyr from Tendulkar et aL| ( 2013 ) 
where the uncertainties are much smaller. 


The X-ray observations with XMM-Newton (Sasaki et al 


|2004| ) show two prominent features of the SNR. Firstly, it has a 
semicircular shape due to a giant molecular cloud (GMC) com¬ 
plex in the west and secondly, a bright diffuse emission region, 
which is known as t he “Lobe”, can be seen (see Fig. (right)). 
It was discovered by |Tatematsu et al.| ( |1987| ) that the GMC com¬ 
plex extends even to the front of the remnant. Consequently, the 
Lobe could simply result from a projection effect, i.e. the result 
of a “hole”, a region with little or no absorption for X-rays, in 
the GMC. Alternatively, this particular emission feature could 
be caused by a shock heated cloud fragment, yielding bright X- 
ray emission. 

In [Sasaki et al.| ( |2006| ) the velocity profiles (from high- 
resolution data from the Five College Radio Astron¬ 

omy Observatory) and Chandra data of three molecular clouds 
around the Lobe were analysed. One of the clouds partly overlap 
with the Lobe. In various parts of this cloud the ^^CO velocity 
profile fits very well with a Gaussian. However, where the Lobe 
and the cloud overlap, the velocity profile deviates from being 
Gaussian and has an additional component towards higher neg¬ 
ative velocities. This indicates an additional acceleration in this 
part of the cloud, e.g. by a supernova (SN) blast wave. Moreover, 
the molecular hydrogen column densities in this region are rel¬ 
atively high, while the foreground absorption in X-rays is lower 
than in other parts. This could also be explained by an interac¬ 
tion between a SN shock wave and a dense cloud. Consequently, 
these results give new evidence for the hypothesis that the Lobe 
is the result of a shocked dense cloud. 

In this work we will verify this hypothesis by 3D hydrody¬ 
namic simulations of the SNR. For this purpose an initial den¬ 
sity distribution, which is based on observations, is sup¬ 
plemented with an additional cloud. By varying the properties of 
this homogeneous cloud, i.e. its position, size, shape, orientation, 
and density, the observed morphology and the X-ray emission of 
CTB 109 can be reproduced very well. 


2. Hydrodynamic simulations 

2.1. Previous modelling efforts 

2.1.1. Analytic models 

The first efforts to model the evolution of CTB 109 were made 
by [Hughes et al.|(|1981|) w ho used the Sedov-Taylor solution 
( [Sedov |T959t [Taylor[[1946[ [1950| ) to estimate the age of the 
remnant. For a homogeneous spherically symmetric medium 
they obtained an age of 17,000 yr from the Einstein data. In 
1991, when the semicircular morphology of the SNR was al¬ 
ready known, [Ni et al.[ ( [1990[ ) used the same method only for 
the eastern part and obtained an age of 13,000 yr. For the west¬ 
ern part they used the analytic solution for the pressure-driven 
snowplow phase (e.g. [McKee & Ostriker][1977[ ) to describe the 
propagation of the blast wave into the GMC. In [Sasaki et al.[ 
( |2004[ ) and [Sasaki et al.[ ( [2013[ ), observational data again only 


from the eastern part obtained with XMM-Newton and Chan¬ 
dra, respectively, were used to estimate the age of CTB 109 with 
the Sedov-Taylor solution. For the XMM-Newton data an age of 
~ 9,400 + 1,000 yr was derived while the Chandra data yielded 
an age of 14,000 + 2,000 yr. The same age estimation was ob¬ 
tained from Suzaku data by [Nakano et al.[ ( 2015[ ) when applying 
the Sedov-Taylor similarity solution. For the alternative assump¬ 
tion that the SNR is in the snowplow phase, they derive 10,000 yr 
as an age. 


2.1.2. Hydrodynamic models 

A first model of CTB 109 which is not based on a similarity solu¬ 
tion is describe d by [Wang et al.[([1992| ). They used the thin shell 
approximation ( [Kompaneets| 1960 ) to describe the hydrodynam¬ 
ics of the semicircular SNR. Therefore, they divided the system 
into a uniform interstellar medium in the east, in which a SN 
explodes, and a uniform dense molecular cloud in the west. The 
interface between the two media is planar and 2 pc away from 
the explosion centre. Around the explosion centre, the domain 
is divided into 100 rings for which the spherically symmetric 
hydrodynamic equations are solved, using a numerical code by 
[Mac Low & McCray| ( [1988| ) and [Chu & Mac Lo^ ( |1990[ ). This 
simulation reproduces the semicircular morphology of CTB 109 
quite well. For an initial energy of £"0 = 3.6 • 10^^ erg, an am¬ 
bient ISM density of ^0 = 0-13 cm“^, and a cloud density of 
^cioud = 36 cm“^ the shape of the simulated remnant matches the 
Einstein data after 13,000 yr. 

In 2012 a sophisticated broadb and study of CTB 109 was 
carried out by [Castro et al.[ ( [2()12[ ). Their spherically symmet¬ 
ric model includes hydrodynamics, efficient cosmic-ray acceler¬ 
ation, non-thermal emission, and a self-consistent calculation of 


the X-r ay t hermal emission as des cribed in detail by [Ellison et al. 
( [2007| ) and [Patnaude et al.[ ( |2009[ ). The evolution of the remnant 
with a power-law ejecta density profile is considered in a circum- 
stellar medium (CSM) with uniform density. The main focus of 
their study is to reproduce the broadband radio. X-ray, and y-ray 
spectrum of CTB 109. For two distinct parameter sets, describ¬ 
ing the hadronic and leptonic components for the y-ray flux rea¬ 
sonable fits were found. However, because the X-ray spectrum 
of the SNR is dominated by thermal emission, a scenario where 
both relativistic leptons and hadrons contribute to the emission 
is more likely. Taking this into account, they found a better fit 
to the observations, which gives an ambient CSM density of 
no = 0.5 cm“^ and an age of 11,000 yr for a distance to CTB 109 
of 2.8 kpc. 


2.1.3. Comparison with our model 

Previous models of CTB 109 concentrated mainly on the esti- 
mation of the age ([Hughes et al.[[198l| [Ni et al.[|1990t [Sasaki 
et al.|200'4|[2013t[Nakano et al.|2Q15| ), only [Wang et al.[ ( |1992 l 


aimed at reproducing of a simple semicircular morphology, and 
[Castro et al.[ ( [2012[ ) focused on the broadband emission. For 
this purpose, simple spherically symmetric hydrodynamic mod¬ 
els are sufficient. However, aiming at a detailed simulation of 
the observed complex morphology, which differs from being 
purely semicircular, no spherically symmetric model can be 
used. Therefore, we use a highly inhomogeneous ambient den¬ 
sity structure from observation, which includes the GMC in the 
west. This would favour a 2D simulation. However, developing 
turbulent structures and the occurring instabilities, which yield 
the break-up of the SN shell, are entirely different in 2D than 
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RA (J2000.0) 


Fig. 1 . Left: 1420-MHz radio continuum image of CTB 109 ( [Kothes et al.|2002 i from the Canadian G alactic Plane Surve y (CGPS, Taylor et al. 
|2003| ). Right: Intensity map (0.3-4.0keV) of CTB 109 in false colour from the XMM-Newton EPIC data ( [Sasaki et al.|2004| ). The very bright point 
source is the AXP IE 2259+586 and the diffuse emission at RA = 23^02"^, Dec = +58°55' (J2000.0) with an extent of ~ 7' is the Lobe. Chandra’s 
field of view is marked as a white rectangle. 


in 3D ( jElmegreen & Scal^|2QQ4t [Kane et al.||200Q| ). While tur¬ 
bulence in nature is an inherently 3D process, one might think 
of turbulent flows that are approximately 2D, in the sense that 
large scale coupling of eddies clearly exceeds the extension in 
one particular dimension. However, the physical properties be¬ 
tween 2D and 3D turbulence are vastly different, so that sacrific¬ 
ing one dimension for a higher spatial resolution might lead to 
erroneous r esult s. In contras t to 3D , as has been noted by [Batch-] 
[elor[ ( [I969[ ) and [Kraichnan] \\961) , due to the conservation of 
vorticity in 2D (in the incompressible, inviscid case without ex¬ 
ternal forcing), the energy cascade is reversed, i.e. from small 
to large scales, while enstrophy flows in the opposite direction 
(double-cascade). In particular, we want to model the bright X- 
ray emission feature inside the SNR for the first time as a result 
of a shocked dense cloud. For this reason an additional cloud is 
added to the initial density structure, precluding again a ID or 
2D model. 

To estimate the age of CTB 109 and to reproduce its semi¬ 
circular morphology, a complex 3D model is not necessary and 
make it even more difficult to study many different parameter 
sets. However, our simulations complement the previous stud¬ 
ies on CTB 109 by analysing the effect of an inhomogeneous 
medium on the evolution of the SNR and by simulating the bright 
X-ray emission region inside the SNR by a shocked dense cloud. 
These aspects could not be studied with their spherically sym¬ 
metric models. 


The details of the thermal and non-thermal emission of the 
SNR, which are modelled in detail by [Castro et~ar] ( [2012[ ), are 
beyond the scope of our study. We investigate in detail the X-ray 
emissivity to compare our finding with the X-ray observations. 
In forthcoming studies, spectral modelling will also be consid¬ 
ered. 


2.2. Numerical details of our studies 

Reproducing observations by numerical simulations represents 
an inverse problem, which, as it is well-known, does not always 
possess a unique solution or even a solution at all. We simplify 


the problem by assuming that the background medium, due to its 
high inertia, has not changed significantly during the short evo¬ 
lution of the SNR, or in other words, we use the present ambient 
medium as the initial model. 


The numerical simulations were performed with the RAM¬ 
SES code ( [Teyssier[2002[ ). The code works on ID, 2D, or 3D un¬ 
structured Cartesian grids and solves the discretised Euler equa¬ 
tions in their conservative form with a second-order Godunov 
scheme for a perfect gas (Piecewise Linear Method). In particu¬ 
lar, the MUSCL scheme (Monotonic Upstream-Centred Scheme 
for Conservatio n Law, [van Leer [1976| \\919) with the HLLC 
Riemann solver ( Toro et al.|1994 ) was used where the MinMod 
slope limiter ( [Roe| 19 81 [ ) is employed to obtain a total variation 
diminishing discretisation scheme ( [Harten|I983[ ). A 58 pc cube, 
which comprises the SNR and the GMC, is initially uniformly 
discretised into 128^ cells, hence the initial resolution of the grid 
is 0.45 pc. Based on a refinement strategy, which refines cells 
adaptively if steep pressure gradients occur, the grid is refined 
during the evolution of the SNR. The finest resolved structures 
due to the adaptive mesh refinement (AMR) have a size of 0.1 pc. 


Given that the age of the SNR is ~ 10"^ yr, we can omit 
its free-expansion phase in the simulation because the duration 
of this phase is only < 1000 yr. The reverse shock must have 
formed and already reached the centre of the SNR heating the 
ejecta. Therefore, the Sedov-Taylor similarity solution can be 
applied. Even though the ejecta might have had inhomogeneities 
like a Si clump as suggested by [Sasaki et al.[ ( [20T^ , such clumps 
should have no significant effect on the overall evolution of the 
forward shock of the SNR expanding into the ISM. Furthermore, 
the end of the free-expansion phase is characterised by the fact 
that the mass of the swept-up interstellar material exceeds the 
mass of the ejecta, i.e. we can neglect the ejecta in our simu¬ 
lations. Consequently, the SN explosion can be modelled by a 
point-explosion, i.e. the pressure po in the 2^ cells at the ex¬ 
plosion centre is given as po = 3 (y - 1) £’o/(47r(Ax)^) with 
Ax = 0.45 pc which yields the expanding remnant. As long as 
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the dynamical time t is smaller than the cooling time 


^cool — 


3kBT 
2nA(T, n) 


( 1 ) 


(e.g. |Spitzer|1978| ) with number density n, temperature T, cool¬ 
ing function A, and Boltzmann constant the expansion is adi¬ 
abatic. Because the effects of cooling and heating cannot be ne¬ 
glected in interaction regions with dense matter, a cooling func¬ 
tion, which was computed with the freely available spectral syn¬ 
thesis code CLOUDY (|Ferland et al.|19981), is employed in the 
RAMSES code. 

The CLOUDY code simulates the conditions within an as- 
trophysical plasma and then predicts the emitted spectrum. In 
particular, the ionisation level in the plasma is determined by 
balancing all ionisation processes (photo. Auger, and colli- 
sional ionisation as well as charge transfer) and all recombi¬ 
nation processes (radiative, low-temperature dielectronic, high- 
temperature dielectronic, three-body recombination, and charge 
transfer) of the 30 lightest elements H - Zn in all stages. For 
this purpose, an optically thick gas slab is divided into a large 
number of zones, such that constant conditions can be assumed 
in each zone. For the free electrons in each zone it is assumed 
that their velocity is predominately Maxwellian distributed with 
a kinetic temperature determined by the balance between heating 
(e.g. photoelectric, mechanical, cosmic-ray) and cooling (e.g. 
forbidden-line cooling and inelastic collisions between electrons 
and other particles). Simultaneously, the associated line and con¬ 
tinuum radiative transfer is solved. The resulting cooling func¬ 
tion for solar abundances, dependent on the number density and 
the temperature, is then used by the RAMSES cooling routine. 


2.3. Model setup 

For a realistic 3D numerical simulation of a supernova ex¬ 
plosion, which reproduces a geometry comparable to that of 
CTB 109, we use the emission of the Canadian Galactic 
Plane Survey (CGPS, [Taylor et al.|2003| ) at a radial velocity of 
V = -51 km s"^ as an indicator of the inhomogeneous ambient 
density structure. This reflects the fact that the molecular mate¬ 
rial as sociated with CTB 109 has approximately this radial ve¬ 
locity ( jTatematsu et al.|1990| ). The idea of using this as a density 
indicator is based on the assumption that the GMC complex in 
the west was hardly displaced by the shock wave of the explo¬ 
sion (initial explosion energy Eq = 10 ^^ erg) because of its high 
inertia. 

For the conversion of the ^^CO data to an ambient density 
structure, results from XMM-Newton (jSasaki et al.||2004] ) and 
Chandra observations ( [Sasaki etal. [20 13 ), respectively, are used. 
For this purpose, the mean ^^CO emission is associated in each 
case with the derived preshock density value, yielding a conver¬ 
sion factor to convert the ^^CO emission to a density structure. 
The XMM-Newton observations aimed at the studies of the mor¬ 
phology CTB 109 and any morphological connection between 
the AXP and the Lobe, which was not found. Moreover, only 
little spectral variation across the remnant were found. The ob¬ 
tained spectra of the shell and the Lobe could be well fitted 
with a one-component non-equilibrium ionisation (NEI) model, 
describing a mixture of the shocked ISM and shocked ejecta. 
The preshock density of the ambient medium was derived to be 
^0 = (0.16 + 0.02) cm“^, where J 3 = D/3.0 kpc is a scal¬ 

ing factor accounting for the distance D of CTB 109, which was 
estimated as 3. 2 + 0.2 kpc by|Kothes & Foster ( 2012| ) with the 
new method of [Foster & Mac Williams | ( [2006| ). This results in a 
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Fig. 2. Density structure of CTB 109, obtained by converting the mean 
^^CO emission in the eastern part with uq - 0.155 cm“^ from XMM- 
Newton data. Overlaid black contour lines display th e 1420-MHz con - 
tinuum emission for 8 K, 15 K, 21 K, 29 K, and 100 K ( jTianet al.|2010j ). 
The current position of the associated AXP IE 2259+586 is marked 
with a magenta star. The circumference of a foreground cloud fragment, 
which was removed for the calculations, is marked as a blue ellipse. 


density of the medium of = 0.155 cm ^ in case of the XMM- 
Newton data. 

The study by [Sasaki et al. ( 2Q13| ) of the Chandra data is com¬ 
plementary to the XMM-Newton study. Only the Lobe and the 
north-eastern part of CTB 109 were observed to study the Lobe 
in greater detail. In comparison to the XMM-Newton observa¬ 
tions, more point sources were found and structures, for exam¬ 
ple in the Lobe, could be resolved. For the spectral analysis the 
emission region was divided into small regions with similar sur¬ 
face brightness resulting in a higher spatial resolution than for 
the XMM-Newton observations. In 50% of the small regions a 
two-component NEI model, describing the shocked ISM and the 
shocked ejecta separately, led to significantly better fits than a 
one-component NEI model. For the density of the ISM a median 
value ^1,median = L2 cm“^ was found. Assuming a compression 
ratio of 4 (y = 5/3), this yields n^ = 0.3 cm“^ as preshock den¬ 
sity in case of the Chandra data. 

We have performed our simulations for both possibilities of 
the ambient medium, i.e. for a mixture of ISM and ejecta as a 
single fluid with n^ = 0.155 cm“^ {XMM-Newton) and only for 
the ISM with n^ = 0.3 cm“^ {Chandra), respectively. As we are 
aiming to explain the Lobe as a shocked dense cloud, the shape, 
size, orientation, and position of this Lobe-forming cloud should 
be independent from the assumption of the mean preshock den¬ 
sity. 

Both density values were used separately to convert the ^^CO 
emission of the CGPS into a density structure by assuming that 
a typical n^ corresponds to the mean value of the ^^CO emission 
in the eastern part (see Fig. I3- 

For the hydrodynamic simulations the 3D initial density dis¬ 
tribution is modelled as a cylinder with the long axis of the 
cylinder running along the line of sight, while its cross section 
is shown in Fig. [^ Initially, for the thermal pressure a uniform 
temperature value of 20 K is adopted and the entire structure is 
considered to be at rest. 

We tested the hypothesis that the Lobe is the result of an 
encounter of a shock wave with a dense cloud. A simply shaped 


Article number, page 4of| 

































J. Bolte et aL: 3D hydrodynamic simulations of the Galactic supernova remnant CTB 109 


Table 1. Variation ranges of the parameter studies to determine the ini¬ 
tial properties of the cloud, which develops into the Lobe. 


Varied quantity 

Range 

Minimum step size 

position i 

109.2°-109.1° 

0.01° 

position h 

(-1.1°)-(-0.9°) 

0.01° 

maximum size 

0.45-8.1 pc 

0.45 pc 

shape 

spherical, elliptical 


orientation 

(-90°) - (-l90°) 

5° 

density 

1-60 cm“^ 

0.5 cm“^ 



109.5° 109.3° 109.1° 108.9° 108.7° 108.5° 


Galactic Longitude 

Fig. 3. Initial density structure from XMM-Newton data, which gives 
the best fit of the 1420-MHz contours of CTB 109 and the Lobe region 
after 8,000 yr. The position of the supernova explosion is marked with a 
black star and the position of the AXP with a magenta star, which is al¬ 
most completely covered by the black star. The elliptically shaped cloud 
(^cioud = 20 cm"^), which resembles the Lobe region after its encounter 
with the shock wave, is located at b) = (109.1545°, -1.0078°). 


the explosion is (£,b) = (109.0962°, -0.9942°), also marked in 

Fig-IU 

For the preshock density value from Chandra observations 
(^0 = 0-3 cm“^) our calculation yields a remnant age of 
11,000 yr and hence the most probable site of the explosion is 
{l,b) = (109.0995°,-0.9936°). 

For both preshock density values the dense homogeneous 
cloud, which becomes the Lobe region, is modelled as an ellip¬ 
tical cylinder of height 1.81 pc at (A b) = (109.1545°, -1.0078°) 
with a semi-major axis of 2.27 pc and a semi-minor axis of 
2.04 pc, which is oriented 5° towards the east. In the case of the 
XMM-Newton data the cloud has a density of ^cioud = 20 cm“^ 
while for the Chandra data the density is ^cioud = 25 cm“^. 
We note that the modelled final structure of the X-ray feature 
depends quite sensitively on the initial properties of the Lobe¬ 
forming cloud, while the morphology of the remnant remains 
almost unaffected. For instance, a higher density leads to a ring¬ 
like emission region, whereas a lower density leads only to a 
diffuse emission. Moreover, the derived age of the SNR also de¬ 
pends on the detailed properties of the cloud; for example a cloud 
which is cut in half, but has the same total mass as before, yields 
an age of the SNR of ~ 10,000 yr for n^ = 0.155 cm“^ because 
the evolution of the feature has to be followed for a longer time 
to give good agreement with the Hi data. 


3.2. Resulting structures of the SNR 

In the following we discuss simulations for both the XMM- 
Newton and Chandra data in order to assess the variations of 
the model parameters compatible with existing observations. 


3.2.1. Estimation of the age 


cloud (resembling the observed Lobe region after its encounter 
with the shock wave) is introduced ad hoc into the initial density 
distribution. 


3. Results 


3.1. Derived initial properties of the Lobe-forming cloud 

Based on the observed shape of the Lobe and the morphol¬ 
ogy of the remnant a good fit for the initial density models 
from XMM-Newton and Chandra data, respectively, could be 
obtained (see Fig. for the model from XMM-Newton data) 
by varying the position, size, shape, orientation, and density of 
a simply shaped cloud. Especially, the position, size, and den¬ 
sity are very sensitive parameters, which significantly influence 
the results. Therefore, we put special emphasis on determin¬ 
ing these quantities. Over 800 different model configurations 
were tested for both mean preshock densities, for details on 
the parameter studies see Table Moreover, the original po¬ 
sition of the supernova explosion is derived from the remnant 
age obtained from the simulation and the proper motion of the 
AXP IE 2259-I-586, which is associated with the SNR CTB 109. 
In [Tendulkar et al.| p013| ) its proper motion was derived as 
(Ma,Ps) = (-6.4 + 0.6 masyr“^ -2.3 + 0.6 masyr“^). 

Eor no = 0.155 cm“^ (XMM-Newton data) our calculation, 
based on the cloud properties shown in Eig. yields a remnant 
age of 8,000 yr, which is resulting in the most probable site of 


Assuming the preshock density distribution derived from the 
XMM-Newton data and a supernova that occurred at (i,b) = 
(109.0962°, -0.9942°), the resulting density structure after 
8,000 yr matches very well the 1420-MHz contours of CTB 109 
in the eastern and southern region as illustrated in Eig. (left). 
Eor the increased preshock density from the Chandra data and 
the site of the explosion at (L,b) = (109.0995°,-0.9936°), the 
resulting density structure after 11,000 yr (see Eig. right) 
matches again very well the 1420-MHz contours of CTB 109 
in the eastern and southern region. 

In summary, we can infer that the evolution time of the SNR 
is in general much smaller than the cooling time for both models, 
such that SNR CTB 109 is still in its Sedov-Taylor phase. How¬ 
ever, this is of course not the case for the Lobe-forming cloud 
and the interaction region with the CMC where cooling cannot 
be neglected. Given the cooling time as in Equation ([T]), we see 
that after ^ ~ 1,000 yr, when the shock has overrun the Lobe¬ 
forming cloud and has heated it up, the ratio ^/^cooi becomes 
greater than 1. The cooling time in the cloud remains thereafter 
at ~ 1,000 yr. The same happens after r ~ 4,000 yr in the west¬ 
ern part where the GMC is. There the cooling time goes even 
down to ~ 100 yr. In all other parts of the SNR ^cooi varies be¬ 
tween ~ 10^ yr and >10^ yr, such that the expansion into the 
eastern ambient medium is adiabatical. We note that simulations 
of the remnant without cooling gave almost identical results for 
CTB 109, which supports our conclusion that the SNR is in its 
Sedov-Taylor phase. 
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Fig. 4. Left: Resulting density distribution from the model fit after 8,000 yr for the XMM-Newton data with the 1420-MHz contours of CTB 109 
overlaid in black. Right: Resulting density distribution from the model fit after 11,000 yr for the Chandra data with the 1420-MHz contours of 
CTB 109 overlaid in black. The position of the supernova explosion is marked with a black star and the position of the AXP with a magenta star. 


3.2.2. Computed emissivities 

To compare the numerical simulations with the X-ray intensity 
from XMM-Newton observations ( [Sasaki et al.||2004| ), we com¬ 
puted the X-ray emissivity £ = A(T, n) from both simulations, 
where A denotes the X-ray cooling function from CLOUDY. 

In Fig. (left) the simulated X-ray intensity for the preshock 
density obtained from XMM-Newton data is displayed. The sim¬ 
ulated X-ray intensity for the Chandra density model is shown 
in Fig.|^(right). 

It can be seen by comparison with Fig. [T] (right) that in 
both models the bright X-ray emission in the Lobe region can 
be reproduced very well, while its morphology is only poorly 
matched by our model assumption of a shocked and simply 
shaped additional cloud. For a better agreement with the ob¬ 
served Lobe structure the additional cloud needs to be inhomo¬ 
geneous and of complex shape resulting in a very large number 
of free parameters, which is beyond the scope of this study. 


3.2.3. Computed velocities 


In Fig. the computed radial velocity from the explosion cen¬ 
tre and the supersonic Mach numbers of our fiducial models are 
displayed. Regions with Mach numbers less than 1 are not dis¬ 
played, hence they appear white. The images in the left column 
show the results for which we assumed the initial conditions de¬ 
rived from the XMM-Newton data, and the right column shows 
the results for the Chandra case. 


The average of the simulated shock velocities for the XMM- 
Newton case (Fig. (upper left) is consiste nt with the es t imated 
mean shock velocity of 720 + 60kms“^ ^ 


by Sasaki et al. 
3 (Fig.|6 


2004| l. 


(lower 


In the region of the shock the Mach number is 
left), which is at the limit of a strong shock. In the region of the 
shocked additional cloud the Mach number exceeds the value of 
5. 


In particular, two regions with high Mach numbers can be 
identified. Firstly, there is a “blob” in the eastern part where the 
shell breaks up and hot gas (~ 10^ K) escapes with high veloc¬ 
ity from the interior of the SNR. Secondly, the Mach number 
is higher in a shell-like structure where the dense Lobe-forming 
cloud slowly dissolves. 
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The computed radial velocity from the explosion centre and 
the supersonic Mach numbers of our fiducial model for the 
Chandra data are displayed in Fig. (right). The mean value 
of the simulated shock velocities is again in agree ment with the 


Sasaki et al. 


estimated mean shock velocity of 460+30 km s“^ by 
(2013). The value of the Mach number is again ~ 3 in the region 
of the shock and also exceeds the value of 5 where the additional 
cloud was shocked. The same high Mach number structures as 
in the case of the XMM-Newton data can also be seen here. 


4. Discussion 


Observation s of the SNR CTB 109, e.g. by [Sasaki et al.| ( [2004 


[2006[ [2013[ ), give evidence that the bright diffuse X-ray emis¬ 


sion in the remnant is the result of a shocked dense cloud. Our 
hydrodynamic simulations support this hypothesis. 

For our best fit model of the initial cloud properties we ob¬ 
tained a density of ^cioud = 20 cm“^ {XMM-Newton data) and 
'^cioud = 25 cm“^ {Chandra data), respectively, and a maximum 
size d = 4.54 pc for the shocked cloud. If we take into account 
that the ambient ISM density derived from the XMM-Newton and 
the Chandra spectra have uncertainties between 10 to 60% be¬ 
cause of uncertainties in the spectral fits and differences in the 
spectra for different extraction regions, the cloud densities cal¬ 
culated for the XMM-Newton and the Chandra case are con¬ 
sistent. The obtained values are much higher than the findings 
of the simple analysis presented in [Sasaki et al.[ ( [2004| ), namely 
^cioud ^ 5 cm“^ and J 1 pc, which was based on the analytic 
estimates given by Sgro[([1975[). How ever, the density values are 
lower than the value of [Castro et ar| ( [2012[ ) from Fermi observa¬ 
tions, which is ^cioud ~ 120 cm“'^. 

Our simulation reproduces for the preshock density n^ = 
0.155 cm“^ (mixture of ISM and ejecta) from XMM-Newton ob¬ 
servations, an initial energy of £”0 = 10^^ erg, and the SN posi¬ 
tion {i,b) = (109.0962°,-0.9942°) after 8,000 yr of evolution 
the observed morphology and size of CTB 109 (cf. Fig.|^. Sim¬ 
ilarly, the morphology of the remnant is matched very well af¬ 
ter 11,000 yr for the higher value n^ = 0.3 cm“^ (ISM only) 
obtained from the Chandra data. These positions are consistent 
with the position of the SN that can be derived from the proper 
motions measurements for the AXP by [Tendulkar et aL (2013). 
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Fig. 5. Left: Modelled X-ray emissivity s - rf' A(T,n) after 8,000 y r (XMM-Newton density model) simulating the intensity map (0.3-4.0keV) 
of CTB 109 from the XMM-Newton EPIC data ( |Sasa^ et al.||2004| (cf. Fig. (right)). Right: Modelled X-ray emissivity s - A(T,n) after 
11,000 yr {Chandra density model) using the same intensity map (0.3-4.0keV) of CTB 109 from the XMM-Newton EPIC data. The position of 
the supernova explosion is marked with a black star and the position of the AXP with a magenta star. 


The SN position, wh ich can be calculated from the proper mo¬ 
tion measurement by [Kaplan et al.| ( |2009| ) is rather inconsistent 
with our calculations. 


The obtained age is comparable with the results from|Sasa^ 


et al.| ( [2004[ ) if the lower preshock density value based on the 


XMM-Newton data is taken. For J 3 = 3.2/3.0 they estimated an 
initial energy of Eq = (8.7 + 3.4) • 10^^ erg and a remnant age of 
t = (9.39+0.96)-10^ yr. An analysis of the C handra observation s 
gives a higher value oft = {\A± 0.2) • 10"^ yr (Sasaki et al. 2013|, 
which was recently verified by Suzaku observations ( [Nakanq 
et al.||2015) . This is comparable with the 11,000 yr, which we 


obtained from our simulations with the Chandra data, and also 
with the age of the remnant, which [Wang et aL] \\992) derived. 
They have modelled CTB 109 numerically with the thin shell ap¬ 
proximation and assumed an initial energy of Eq = 3.6 • 10^^ erg, 
an ambient ISM density of ^0 = 0-13 cm“^, and a cloud density 
of ^cioud = 36 cm“^. After t = 1.3 • 10"^ yr their model repro¬ 
duced a semicircular shell of the observed size. In contrast to 
their thin shell approximation, our model does not need a spe¬ 
cial geometry or symmetry assumptions, but solves the hydro- 
dynamic equations for a completely inhomogeneous, complex 
3D realistic ISM with high resolution. Consequently, our model 
gives improved estimates on the initial energy and the remnant 
age and can also explain its observed morphology. 


In a similar manner we see that the shock velocities from the 
simulations and observations also agree very well for the differ¬ 
ent preshock densities. For the simulation with the XMM-Newton 
data the mean shock velocity fits with = 720 ± 60kms“^ 
( [Sasaki et al.||2004| ) and the simulations wi th the Chandra dat a 
agree very well with = 460 ± 30kms“^ ( Sasaki et al.[|2013' 


The latter value of shock velocity was also obtained from the 
recent Suzaku observations ( [Nakano et al.|2015] ). 


A high synchrotron emission, which is visible in the ra¬ 
dio continuum, correlates with more efficiently accelerated elec¬ 
trons. Hence, higher shock velocities are expected in this region. 
The distribution of the shock velocities in Fig. show high ve¬ 
locities in a shell-like structure in the north-east and more con¬ 
centrated in the south. These structures match quite well with 
the 1420-MHz radio continuum of CTB 109, displayed in Fig.[2 


(left). This indicates that the simulated hydrodynamic structure 
is a plausible representation of the remnant’s structure. 

At present, all parameter studies that we have performed 
give very good agreement with both XMM-Newton and Chan¬ 
dra data sets. Therefore, we cannot favour one model over the 
other. However, since the Chandra data have been fitted with a 
more realistic two-component NEI model instead of a simplistic 
one-component NEI model, we believe that the density model 
obtained from the Chandra data is a better description of the re¬ 
ality. Consequently, the obtained cloud properties, remnant age, 
and shock velocities from this model should be more reliable. 


5. Concluding remarks 


The presented numerical simulations of a supernova explosion 
in a realistic inhomogeneous medium make it very plausible that 
the Lobe region is the result of a shocked dense elliptical cloud. 
In the best fit initial density model the cloud has a density of 
^cioud = 25 cm“^, which lies between the values from the simple 
analysis by [Sasaki et al.[ ( [2004| ) and the estimate by [Castro et al.[ 
( [2012| ). Not only the observed morphology of the remnant and of 
the Lobe, but also the distribution of the observed X-ray intensity 
in the Lobe region could be reproduced very well by this model. 
For an initial explosion energy of £”0 = 10^^ erg the remnant 
age can be estimated to be 11,000 yr. The result agrees with the 
remnant age of ~ 14,000 yr derived by [Sasaki et al.[ ( [2013| ). 
This also allows us to derive the most probable site of the SN 
explosion. 

Our simulations have shown that present observations in 
^^CO, Hi, and X-ray, constrain the initial conditions well enough 
to derive important physical parameters, e.g. explosion energy 
and also evolution time of an SNR. The simulations of the SNR 
CTB 109 show also that the preshock density distribution must 
have been more complex than simply a medium with a homoge¬ 
neous density gradient. 

Thus, we are confident that our method of using the present 
ambient density distribution of the SNR and modifying it gradu¬ 
ally to reproduce the morphology and X-ray emissivity of young 
SNRs by an ab initio 3D hydrodynamic simulation is reasonable 
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Fig. 6. Left: Radial velocity from the explosion centre after 8,000 yr of our fiducial model based on the XMM-Newton data (upper panel). Supersonic 
Mach numbers after 8,000 yr for our fiducial model based on the XMM-Newton data (lower panel). Right: Radial velocity from the explosion centre 
after 11,000 yr of our fiducial model based the Chandra data (upper panel). Supersonic Mach numbers after 11,000 yr for our fiducial model based 
on the Chandra data (lower panel). Regions with Mach numbers less than 1 are not displayed, hence they appear white. The position of the 
supernova explosion is marked with a black star and the position of the AXP with a magenta star. 


and will help to understand the realistic conditions. This method 
will be explored in detail in forthcoming studies. 
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